Let us set some global options for all code chunks in this document.

# Set seed for reproducibility
set.seed(1982) 
# Set global options for all code chunks
knitr::opts_chunk$set(
  # Disable messages printed by R code chunks
  message = FALSE,    
  # Disable warnings printed by R code chunks
  warning = FALSE,    
  # Show R code within code chunks in output
  echo = TRUE,        
  # Include both R code and its results in output
  include = TRUE,     
  # Evaluate R code chunks
  eval = TRUE,       
  # Enable caching of R code chunks for faster rendering
  cache = FALSE,      
  # Align figures in the center of the output
  fig.align = "center",
  # Enable retina display for high-resolution figures
  retina = 2,
  # Show errors in the output instead of stopping rendering
  error = TRUE,
  # Do not collapse code and output into a single block
  collapse = FALSE
)
# inla.upgrade(testing = TRUE)
# remotes::install_github("inlabru-org/inlabru", ref = "devel")
# remotes::install_github("davidbolin/rspde", ref = "devel")
# remotes::install_github("davidbolin/metricgraph", ref = "devel")
library(INLA)
library(inlabru)
library(rSPDE)
library(MetricGraph)
library(grateful)
library(ggplot2)
library(reshape2)
library(plotly)

We want to solve the fractional diffusion equation \[\begin{equation} \label{eq:maineq} \partial_t u+(\kappa^2-\Delta_\Gamma)^{\frac{\alpha}{2}} u=f \text { on } \Gamma \times(0, T), \quad u(0)=u_0 \text { on } \Gamma, \end{equation}\] where \(u\) satisfies the Kirchhoff vertex conditions \[\begin{equation} \label{eq:Kcond} \left\{\phi\in C(\Gamma)\;\Big|\; \forall v\in V: \sum_{e\in\mathcal{E}_v}\partial_e \phi(v)=0 \right\} \end{equation}\]

If \(f=0\), then the solution is given by \[\begin{equation} \label{eq:sol_reprentation} u(s,t) = \displaystyle\sum_{j\in\mathbb{N}}e^{-\EIGENVAL^{\frac{\alpha}{2}}_jt}\left(u_0, e_j\right)_{L_2(\Gamma)}e_j(s). \end{equation}\]

# Function to build a tadpole graph and create a mesh
gets_graph_tadpole <- function(h){
  edge1 <- rbind(c(0,0),c(1,0))
  theta <- seq(from=-pi,to=pi,length.out = 10000)
  edge2 <- cbind(1+1/pi+cos(theta)/pi,sin(theta)/pi)
  edges = list(edge1, edge2)
  graph <- metric_graph$new(edges = edges)
  graph$build_mesh(h = h)
  return(graph)
}

Let \(\Gamma_T = (\mathcal{V},\mathcal{E})\) characterize the tadpole graph with \(\mathcal{V}= \{v_1,v_2\}\) and \(\mathcal{E}= \{e_1,e_2\}\) as specified in Figure \(\ref{Interval.Circle.Tadpole}\)c. The left edge \(e_1\) has length 1 and the circular edge \(e_2\) has length 2. As discussed in Subsection \(\ref{subsec:prelim}\), a point on \(e_1\) is parameterized via \(s=\left(e_1, t\right)\) for \(t \in[0,1]\) and a point on \(e_2\) via \(s=\left(e_2, t\right)\) for \(t\in[0,2]\). One can verify that \(-\Delta_\Gamma\) has eigenvalues \(0,\left\{(i \pi / 2)^2\right\}_{i \in \mathbb{N}}\) and \(\left\{(i \pi / 2)^2\right\}_{2 i \in \mathbb{N}}\) with corresponding eigenfunctions \(\phi_0\), \(\left\{\phi_i\right\}_{i \in \mathbb{N}}\), and \(\left\{\psi_i\right\}_{2 i \in \mathbb{N}}\) given by \(\phi_0(s)=1 / \sqrt{3}\) and \[\begin{equation*} \phi_i(s)=C_{\phi, i}\begin{cases} -2 \sin (\frac{i\pi}{2}) \cos (\frac{i \pi t}{2}), & s \in e_1, \\ \sin (i \pi t / 2), & s \in e_2, \end{cases}, \quad \psi_i(s)=\frac{\sqrt{3}}{\sqrt{2}} \begin{cases} (-1)^{i / 2} \cos (\frac{i \pi t}{2}), & s \in e_1, \\ \cos (\frac{i \pi t}{2}), & s \in e_2, \end{cases}, \end{equation*}\] where \(C_{\phi, i}=1\) if \(i\) is even and \(C_{\phi, i}=1 / \sqrt{3}\) otherwise. Moreover, these functions form an orthonormal basis for \(L_2(\Gamma_T)\).

# Function to compute the eigenfunctions 
tadpole.eig <- function(k,graph){
x1 <- c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[graph$mesh$PtE[,1]==1,2]) 
x2 <- c(0,graph$get_edge_lengths()[2]*graph$mesh$PtE[graph$mesh$PtE[,1]==2,2]) 

if(k==0){ 
  f.e1 <- rep(1,length(x1)) 
  f.e2 <- rep(1,length(x2)) 
  f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
  f = list(phi=f1/sqrt(3)) 
  
} else {
  f.e1 <- -2*sin(pi*k*1/2)*cos(pi*k*x1/2) 
  f.e2 <- sin(pi*k*x2/2)                  
  
  f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
  
  if((k %% 2)==1){ 
    f = list(phi=f1/sqrt(3)) 
  } else { 
    f.e1 <- (-1)^{k/2}*cos(pi*k*x1/2)
    f.e2 <- cos(pi*k*x2/2)
    f2 = c(f.e1[1],f.e2[1],f.e1[-1],f.e2[-1]) 
    f <- list(phi=f1,psi=f2/sqrt(3/2))
  }
}

return(f)
}

Implementation of \(u\)

h <- 0.001
graph <- gets_graph_tadpole(h = h)
# Compute the FEM matrices
graph$compute_fem()
x <- graph$mesh$V[, 1]
y <- graph$mesh$V[, 2]
edge_number <- graph$mesh$VtE[, 1]
pos <- sum(edge_number == 1)+1
order_to_plot <- function(v)return(c(v[1], v[3:pos], v[2], v[(pos+1):length(v)], v[2]))
weights <- graph$mesh$weights
kappa <- 1
n_finite <- 100

EIGENVAL <- c()       # initialize empty vector for eigenvalues
EIGENFUN <- NULL       # initialize NULL for eigenfunctions matrix
INDEX <- c()

for (j in 0:n_finite) {
    lambda_j <- kappa^2 + (j*pi/2)^2
    e_j <- tadpole.eig(j,graph)$phi#/sqrt(1000)
    EIGENVAL <- c(EIGENVAL, lambda_j)         # append scalar to vector
    EIGENFUN <- cbind(EIGENFUN, e_j)            # append column to matrix
    INDEX <- c(INDEX, j)
    if (j>0 && (j %% 2 == 0)) {
      lambda_j <- kappa^2 + (j*pi/2)^2
      e_j <- tadpole.eig(j,graph)$psi#/sqrt(1000)
      EIGENVAL <- c(EIGENVAL, lambda_j)         # append scalar to vector
      EIGENFUN <- cbind(EIGENFUN, e_j)            # append column to matrix
      INDEX <- c(INDEX, j+0.1)
    }
}

# Compute inner product matrix
INNER_PRODUCTS <- t(EIGENFUN) %*% (EIGENFUN * weights)

# Assuming INNER_PRODUCTS is a numeric matrix
plot_ly(
  z = INNER_PRODUCTS,
  type = "heatmap",
  colorscale = "Viridis"  # You can try "Hot", "Jet", etc.
) %>%
  layout(
    title = "Inner Product Matrix of Eigenfunctions",
    xaxis = list(title = "Index"),
    yaxis = list(title = "Index")
  )
EIGENFUN <- apply(EIGENFUN, 2, order_to_plot)
x <- order_to_plot(x)
y <- order_to_plot(y)

# Create interactive plot
fig <- plot_ly(x = ~INDEX, y = ~EIGENVAL, type = 'scatter', mode = 'markers',
               marker = list(size = 6),
               name = "Eigenvalues")

fig <- fig %>% layout(title = "Eigenvalues",
                      xaxis = list(title = "index i"),
                      yaxis = list(title = "Eigenvalues_i"))



plot_data <- data.frame(
  x = rep(x, times = ncol(EIGENFUN)),
  y = rep(y, times = ncol(EIGENFUN)),
  z = as.vector(EIGENFUN),
  frame = rep(INDEX, each = length(x))
)

# Compute axis limits
x_range <- range(x)
y_range <- range(y)
z_range <- range(EIGENFUN)

# Initial plot setup (first frame only)
p <- plot_ly(plot_data, frame = ~frame) %>%
  add_trace(
    x = ~x, y = ~y, z = ~z,
    type = "scatter3d", mode = "lines",
    name = paste0("n_finite = ", n_finite),
    line = list(color = "blue", width = 2)
  ) %>%
  layout(
    scene = list(
      xaxis = list(title = "x", range = x_range),
      yaxis = list(title = "y", range = y_range),
      zaxis = list(title = "Value", range = z_range)
    ),
    updatemenus = list(
      list(
        type = "buttons", showactive = FALSE,
        buttons = list(
          list(label = "Play", method = "animate",
               args = list(NULL, list(frame = list(duration = 100, redraw = TRUE), fromcurrent = TRUE))),
          list(label = "Pause", method = "animate",
               args = list(NULL, list(mode = "immediate", frame = list(duration = 0), redraw = FALSE)))
        )
      )
    ),
    title = "Time: 0"
  )

# Convert to plotly object with frame info
pb <- plotly_build(p)

# Inject custom titles into each frame
for (i in seq_along(pb$x$frames)) {
  t <- INDEX[i]
  err <- signif(EIGENVAL[i], 4)
  pb$x$frames[[i]]$layout <- list(title = paste0("i: ", t, " | Eigenvalue_i: ", err))
}


fig  # Display the plot
pb
LS0tCnRpdGxlOiAiU29sdmluZyBhIHBhcmFib2xpYyBlcXVhdGlvbiIKZGF0ZTogIkNyZWF0ZWQ6IDIwLTA0LTIwMjUuIExhc3QgbW9kaWZpZWQ6IGByIGZvcm1hdChTeXMudGltZSgpLCAnJWQtJW0tJVkuJylgIgpvdXRwdXQ6CiAgaHRtbF9kb2N1bWVudDoKICAgIG1hdGhqYXg6ICJodHRwczovL2Nkbi5qc2RlbGl2ci5uZXQvbnBtL21hdGhqYXhAMy9lczUvdGV4LW1tbC1jaHRtbC5qcyIKICAgIGhpZ2hsaWdodDogcHlnbWVudHMKICAgIHRoZW1lOiBmbGF0bHkKICAgIGNvZGVfZm9sZGluZzogc2hvdyAjIGNsYXNzLnNvdXJjZSA9ICJmb2xkLWhpZGUiIHRvIGhpZGUgY29kZSBhbmQgYWRkIGEgYnV0dG9uIHRvIHNob3cgaXQKICAgIGRmX3ByaW50OiBwYWdlZAogICAgdG9jOiB0cnVlCiAgICB0b2NfZmxvYXQ6CiAgICAgIGNvbGxhcHNlZDogdHJ1ZQogICAgICBzbW9vdGhfc2Nyb2xsOiB0cnVlCiAgICBudW1iZXJfc2VjdGlvbnM6IGZhbHNlCiAgICBmaWdfY2FwdGlvbjogdHJ1ZQogICAgY29kZV9kb3dubG9hZDogdHJ1ZQphbHdheXNfYWxsb3dfaHRtbDogdHJ1ZQpiaWJsaW9ncmFwaHk6IAogIC0gcmVmZXJlbmNlcy5iaWIKICAtIGdyYXRlZnVsLXJlZnMuYmliCmhlYWRlci1pbmNsdWRlczoKICAtIFxuZXdjb21tYW5ke1xhcn17XG1hdGhiYntSfX0KICAtIFxuZXdjb21tYW5ke1xsbGF2fVsxXXtcbGVmdFx7IzFccmlnaHRcfX0KICAtIFxuZXdjb21tYW5ke1xwYXJlfVsxXXtcbGVmdCgjMVxyaWdodCl9CiAgLSBcbmV3Y29tbWFuZHtcTmNhbH17XG1hdGhjYWx7Tn19CiAgLSBcbmV3Y29tbWFuZHtcVmNhbH17XG1hdGhjYWx7Vn19CiAgLSBcbmV3Y29tbWFuZHtcRWNhbH17XG1hdGhjYWx7RX19CiAgLSBcbmV3Y29tbWFuZHtcV2NhbH17XG1hdGhjYWx7V319Ci0tLQoKYGBge3IgeGFyaW5nYW5FeHRyYS1jbGlwYm9hcmQsIGVjaG8gPSBGQUxTRX0KaHRtbHRvb2xzOjp0YWdMaXN0KAogIHhhcmluZ2FuRXh0cmE6OnVzZV9jbGlwYm9hcmQoCiAgICBidXR0b25fdGV4dCA9ICI8aSBjbGFzcz1cImZhLXNvbGlkIGZhLWNsaXBib2FyZFwiIHN0eWxlPVwiY29sb3I6ICMwMDAwOEJcIj48L2k+IiwKICAgIHN1Y2Nlc3NfdGV4dCA9ICI8aSBjbGFzcz1cImZhIGZhLWNoZWNrXCIgc3R5bGU9XCJjb2xvcjogIzkwQkU2RFwiPjwvaT4iLAogICAgZXJyb3JfdGV4dCA9ICI8aSBjbGFzcz1cImZhIGZhLXRpbWVzLWNpcmNsZVwiIHN0eWxlPVwiY29sb3I6ICNGOTQxNDRcIj48L2k+IgogICksCiAgcm1hcmtkb3duOjpodG1sX2RlcGVuZGVuY3lfZm9udF9hd2Vzb21lKCkKKQpgYGAKCgpgYGB7Y3NzLCBlY2hvID0gRkFMU0V9CmJvZHkgLm1haW4tY29udGFpbmVyIHsKICBtYXgtd2lkdGg6IDEwMCUgIWltcG9ydGFudDsKICB3aWR0aDogMTAwJSAhaW1wb3J0YW50Owp9CmJvZHkgewogIG1heC13aWR0aDogMTAwJSAhaW1wb3J0YW50Owp9Cgpib2R5LCB0ZCB7CiAgIGZvbnQtc2l6ZTogMTZweDsKfQpjb2RlLnJ7CiAgZm9udC1zaXplOiAxNHB4Owp9CnByZSB7CiAgZm9udC1zaXplOiAxNHB4Cn0KLmN1c3RvbS1ib3ggewogIGJhY2tncm91bmQtY29sb3I6ICNmNWY3ZmE7IC8qIExpZ2h0IGdyZXktYmx1ZSBiYWNrZ3JvdW5kICovCiAgYm9yZGVyLWNvbG9yOiAjZTFlOGVkOyAvKiBMaWdodCBib3JkZXIgY29sb3IgKi8KICBjb2xvcjogIzJjM2U1MDsgLyogRGFyayB0ZXh0IGNvbG9yICovCiAgcGFkZGluZzogMTVweDsgLyogUGFkZGluZyBpbnNpZGUgdGhlIGJveCAqLwogIGJvcmRlci1yYWRpdXM6IDVweDsgLyogUm91bmRlZCBjb3JuZXJzICovCiAgbWFyZ2luLWJvdHRvbTogMjBweDsgLyogU3BhY2luZyBiZWxvdyB0aGUgYm94ICovCn0KLmNhcHRpb24gewogIG1hcmdpbjogYXV0bzsKICB0ZXh0LWFsaWduOiBjZW50ZXI7CiAgbWFyZ2luLWJvdHRvbTogMjBweDsgLyogU3BhY2luZyBiZWxvdyB0aGUgYm94ICovCn0KYGBgCgoKTGV0IHVzIHNldCBzb21lIGdsb2JhbCBvcHRpb25zIGZvciBhbGwgY29kZSBjaHVua3MgaW4gdGhpcyBkb2N1bWVudC4KCgpgYGB7cn0KIyBTZXQgc2VlZCBmb3IgcmVwcm9kdWNpYmlsaXR5CnNldC5zZWVkKDE5ODIpIAojIFNldCBnbG9iYWwgb3B0aW9ucyBmb3IgYWxsIGNvZGUgY2h1bmtzCmtuaXRyOjpvcHRzX2NodW5rJHNldCgKICAjIERpc2FibGUgbWVzc2FnZXMgcHJpbnRlZCBieSBSIGNvZGUgY2h1bmtzCiAgbWVzc2FnZSA9IEZBTFNFLCAgICAKICAjIERpc2FibGUgd2FybmluZ3MgcHJpbnRlZCBieSBSIGNvZGUgY2h1bmtzCiAgd2FybmluZyA9IEZBTFNFLCAgICAKICAjIFNob3cgUiBjb2RlIHdpdGhpbiBjb2RlIGNodW5rcyBpbiBvdXRwdXQKICBlY2hvID0gVFJVRSwgICAgICAgIAogICMgSW5jbHVkZSBib3RoIFIgY29kZSBhbmQgaXRzIHJlc3VsdHMgaW4gb3V0cHV0CiAgaW5jbHVkZSA9IFRSVUUsICAgICAKICAjIEV2YWx1YXRlIFIgY29kZSBjaHVua3MKICBldmFsID0gVFJVRSwgICAgICAgCiAgIyBFbmFibGUgY2FjaGluZyBvZiBSIGNvZGUgY2h1bmtzIGZvciBmYXN0ZXIgcmVuZGVyaW5nCiAgY2FjaGUgPSBGQUxTRSwgICAgICAKICAjIEFsaWduIGZpZ3VyZXMgaW4gdGhlIGNlbnRlciBvZiB0aGUgb3V0cHV0CiAgZmlnLmFsaWduID0gImNlbnRlciIsCiAgIyBFbmFibGUgcmV0aW5hIGRpc3BsYXkgZm9yIGhpZ2gtcmVzb2x1dGlvbiBmaWd1cmVzCiAgcmV0aW5hID0gMiwKICAjIFNob3cgZXJyb3JzIGluIHRoZSBvdXRwdXQgaW5zdGVhZCBvZiBzdG9wcGluZyByZW5kZXJpbmcKICBlcnJvciA9IFRSVUUsCiAgIyBEbyBub3QgY29sbGFwc2UgY29kZSBhbmQgb3V0cHV0IGludG8gYSBzaW5nbGUgYmxvY2sKICBjb2xsYXBzZSA9IEZBTFNFCikKYGBgCgoKCgpgYGB7cn0KIyBpbmxhLnVwZ3JhZGUodGVzdGluZyA9IFRSVUUpCiMgcmVtb3Rlczo6aW5zdGFsbF9naXRodWIoImlubGFicnUtb3JnL2lubGFicnUiLCByZWYgPSAiZGV2ZWwiKQojIHJlbW90ZXM6Omluc3RhbGxfZ2l0aHViKCJkYXZpZGJvbGluL3JzcGRlIiwgcmVmID0gImRldmVsIikKIyByZW1vdGVzOjppbnN0YWxsX2dpdGh1YigiZGF2aWRib2xpbi9tZXRyaWNncmFwaCIsIHJlZiA9ICJkZXZlbCIpCmxpYnJhcnkoSU5MQSkKbGlicmFyeShpbmxhYnJ1KQpsaWJyYXJ5KHJTUERFKQpsaWJyYXJ5KE1ldHJpY0dyYXBoKQpsaWJyYXJ5KGdyYXRlZnVsKQpsaWJyYXJ5KGdncGxvdDIpCmxpYnJhcnkocmVzaGFwZTIpCmxpYnJhcnkocGxvdGx5KQpgYGAKCgpXZSB3YW50IHRvIHNvbHZlIHRoZSBmcmFjdGlvbmFsIGRpZmZ1c2lvbiBlcXVhdGlvbgpcYmVnaW57ZXF1YXRpb259ClxsYWJlbHtlcTptYWluZXF9CiAgICBccGFydGlhbF90IHUrKFxrYXBwYV4yLVxEZWx0YV9cR2FtbWEpXntcZnJhY3tcYWxwaGF9ezJ9fSB1PWYgXHRleHQgeyBvbiB9IFxHYW1tYSBcdGltZXMoMCwgVCksIFxxdWFkIHUoMCk9dV8wIFx0ZXh0IHsgb24gfSBcR2FtbWEsClxlbmR7ZXF1YXRpb259CndoZXJlICR1JCBzYXRpc2ZpZXMgdGhlIEtpcmNoaG9mZiB2ZXJ0ZXggY29uZGl0aW9ucwpcYmVnaW57ZXF1YXRpb259ClxsYWJlbHtlcTpLY29uZH0KICAgIFxsZWZ0XHtccGhpXGluIEMoXEdhbW1hKVw7XEJpZ3xcOyBcZm9yYWxsIHZcaW4gVjogXHN1bV97ZVxpblxtYXRoY2Fse0V9X3Z9XHBhcnRpYWxfZSBccGhpKHYpPTAgXHJpZ2h0XH0KXGVuZHtlcXVhdGlvbn0KCklmICRmPTAkLCB0aGVuIHRoZSBzb2x1dGlvbiBpcyBnaXZlbiBieQpcYmVnaW57ZXF1YXRpb259ClxsYWJlbHtlcTpzb2xfcmVwcmVudGF0aW9ufQogICAgICAgIHUocyx0KSA9IFxkaXNwbGF5c3R5bGVcc3VtX3tqXGluXG1hdGhiYntOfX1lXnstXEVJR0VOVkFMXntcZnJhY3tcYWxwaGF9ezJ9fV9qdH1cbGVmdCh1XzAsIGVfalxyaWdodClfe0xfMihcR2FtbWEpfWVfaihzKS4KXGVuZHtlcXVhdGlvbn0KCmBgYHtyfQojIEZ1bmN0aW9uIHRvIGJ1aWxkIGEgdGFkcG9sZSBncmFwaCBhbmQgY3JlYXRlIGEgbWVzaApnZXRzX2dyYXBoX3RhZHBvbGUgPC0gZnVuY3Rpb24oaCl7CiAgZWRnZTEgPC0gcmJpbmQoYygwLDApLGMoMSwwKSkKICB0aGV0YSA8LSBzZXEoZnJvbT0tcGksdG89cGksbGVuZ3RoLm91dCA9IDEwMDAwKQogIGVkZ2UyIDwtIGNiaW5kKDErMS9waStjb3ModGhldGEpL3BpLHNpbih0aGV0YSkvcGkpCiAgZWRnZXMgPSBsaXN0KGVkZ2UxLCBlZGdlMikKICBncmFwaCA8LSBtZXRyaWNfZ3JhcGgkbmV3KGVkZ2VzID0gZWRnZXMpCiAgZ3JhcGgkYnVpbGRfbWVzaChoID0gaCkKICByZXR1cm4oZ3JhcGgpCn0KYGBgCgpMZXQgJFxHYW1tYV9UID0gKFxWY2FsLFxFY2FsKSQgY2hhcmFjdGVyaXplIHRoZSB0YWRwb2xlIGdyYXBoIHdpdGggJFxWY2FsID0gXHt2XzEsdl8yXH0kIGFuZCAkXEVjYWwgPSBce2VfMSxlXzJcfSQgYXMgc3BlY2lmaWVkIGluIEZpZ3VyZSBccmVme0ludGVydmFsLkNpcmNsZS5UYWRwb2xlfWMuIFRoZSBsZWZ0IGVkZ2UgJGVfMSQgaGFzIGxlbmd0aCAxIGFuZCB0aGUgY2lyY3VsYXIgZWRnZSAkZV8yJCBoYXMgbGVuZ3RoIDIuIEFzIGRpc2N1c3NlZCBpbiBTdWJzZWN0aW9uIFxyZWZ7c3Vic2VjOnByZWxpbX0sIGEgcG9pbnQgb24gJGVfMSQgaXMgcGFyYW1ldGVyaXplZCB2aWEgJHM9XGxlZnQoZV8xLCB0XHJpZ2h0KSQgZm9yICR0IFxpblswLDFdJCBhbmQgYSBwb2ludCBvbiAkZV8yJCB2aWEgJHM9XGxlZnQoZV8yLCB0XHJpZ2h0KSQgZm9yICR0XGluWzAsMl0kLiBPbmUgY2FuIHZlcmlmeSB0aGF0ICQtXERlbHRhX1xHYW1tYSQgaGFzIGVpZ2VudmFsdWVzICQwLFxsZWZ0XHsoaSBccGkgLyAyKV4yXHJpZ2h0XH1fe2kgXGluIFxtYXRoYmJ7Tn19JCBhbmQgJFxsZWZ0XHsoaSBccGkgLyAyKV4yXHJpZ2h0XH1fezIgaSBcaW4gXG1hdGhiYntOfX0kIHdpdGggY29ycmVzcG9uZGluZyBlaWdlbmZ1bmN0aW9ucyAkXHBoaV8wJCwgJFxsZWZ0XHtccGhpX2lccmlnaHRcfV97aSBcaW4gXG1hdGhiYntOfX0kLCBhbmQgJFxsZWZ0XHtccHNpX2lccmlnaHRcfV97MiBpIFxpbiBcbWF0aGJie059fSQgZ2l2ZW4gYnkgJFxwaGlfMChzKT0xIC8gXHNxcnR7M30kIGFuZCAKXGJlZ2lue2VxdWF0aW9uKn0KICAgIFxwaGlfaShzKT1DX3tccGhpLCBpfVxiZWdpbntjYXNlc30KICAgICAgICAtMiBcc2luIChcZnJhY3tpXHBpfXsyfSkgXGNvcyAoXGZyYWN7aSBccGkgdH17Mn0pLCAmIHMgXGluIGVfMSwgXFwKXHNpbiAoaSBccGkgdCAvIDIpLCAmIHMgXGluIGVfMiwKICAgIFxlbmR7Y2FzZXN9LApccXVhZCAKICAgIFxwc2lfaShzKT1cZnJhY3tcc3FydHszfX17XHNxcnR7Mn19IFxiZWdpbntjYXNlc30KICAgICgtMSlee2kgLyAyfSBcY29zIChcZnJhY3tpIFxwaSB0fXsyfSksICYgcyBcaW4gZV8xLCBcXApcY29zIChcZnJhY3tpIFxwaSB0fXsyfSksICYgcyBcaW4gZV8yLApcZW5ke2Nhc2VzfSwKXGVuZHtlcXVhdGlvbip9CndoZXJlICRDX3tccGhpLCBpfT0xJCBpZiAkaSQgaXMgZXZlbiBhbmQgJENfe1xwaGksIGl9PTEgLyBcc3FydHszfSQgb3RoZXJ3aXNlLiBNb3Jlb3ZlciwgdGhlc2UgZnVuY3Rpb25zIGZvcm0gYW4gb3J0aG9ub3JtYWwgYmFzaXMgZm9yICRMXzIoXEdhbW1hX1QpJC4KCmBgYHtyfQojIEZ1bmN0aW9uIHRvIGNvbXB1dGUgdGhlIGVpZ2VuZnVuY3Rpb25zIAp0YWRwb2xlLmVpZyA8LSBmdW5jdGlvbihrLGdyYXBoKXsKeDEgPC0gYygwLGdyYXBoJGdldF9lZGdlX2xlbmd0aHMoKVsxXSpncmFwaCRtZXNoJFB0RVtncmFwaCRtZXNoJFB0RVssMV09PTEsMl0pIAp4MiA8LSBjKDAsZ3JhcGgkZ2V0X2VkZ2VfbGVuZ3RocygpWzJdKmdyYXBoJG1lc2gkUHRFW2dyYXBoJG1lc2gkUHRFWywxXT09MiwyXSkgCgppZihrPT0wKXsgCiAgZi5lMSA8LSByZXAoMSxsZW5ndGgoeDEpKSAKICBmLmUyIDwtIHJlcCgxLGxlbmd0aCh4MikpIAogIGYxID0gYyhmLmUxWzFdLGYuZTJbMV0sZi5lMVstMV0sIGYuZTJbLTFdKSAKICBmID0gbGlzdChwaGk9ZjEvc3FydCgzKSkgCiAgCn0gZWxzZSB7CiAgZi5lMSA8LSAtMipzaW4ocGkqayoxLzIpKmNvcyhwaSprKngxLzIpIAogIGYuZTIgPC0gc2luKHBpKmsqeDIvMikgICAgICAgICAgICAgICAgICAKICAKICBmMSA9IGMoZi5lMVsxXSxmLmUyWzFdLGYuZTFbLTFdLCBmLmUyWy0xXSkgCiAgCiAgaWYoKGsgJSUgMik9PTEpeyAKICAgIGYgPSBsaXN0KHBoaT1mMS9zcXJ0KDMpKSAKICB9IGVsc2UgeyAKICAgIGYuZTEgPC0gKC0xKV57ay8yfSpjb3MocGkqayp4MS8yKQogICAgZi5lMiA8LSBjb3MocGkqayp4Mi8yKQogICAgZjIgPSBjKGYuZTFbMV0sZi5lMlsxXSxmLmUxWy0xXSxmLmUyWy0xXSkgCiAgICBmIDwtIGxpc3QocGhpPWYxLHBzaT1mMi9zcXJ0KDMvMikpCiAgfQp9CgpyZXR1cm4oZikKfQpgYGAKCkltcGxlbWVudGF0aW9uIG9mICR1JAoKYGBge3J9CmggPC0gMC4wMDEKZ3JhcGggPC0gZ2V0c19ncmFwaF90YWRwb2xlKGggPSBoKQojIENvbXB1dGUgdGhlIEZFTSBtYXRyaWNlcwpncmFwaCRjb21wdXRlX2ZlbSgpCnggPC0gZ3JhcGgkbWVzaCRWWywgMV0KeSA8LSBncmFwaCRtZXNoJFZbLCAyXQplZGdlX251bWJlciA8LSBncmFwaCRtZXNoJFZ0RVssIDFdCnBvcyA8LSBzdW0oZWRnZV9udW1iZXIgPT0gMSkrMQpvcmRlcl90b19wbG90IDwtIGZ1bmN0aW9uKHYpcmV0dXJuKGModlsxXSwgdlszOnBvc10sIHZbMl0sIHZbKHBvcysxKTpsZW5ndGgodildLCB2WzJdKSkKd2VpZ2h0cyA8LSBncmFwaCRtZXNoJHdlaWdodHMKYGBgCgoKCmBgYHtyfQprYXBwYSA8LSAxCm5fZmluaXRlIDwtIDEwMAoKRUlHRU5WQUwgPC0gYygpICAgICAgICMgaW5pdGlhbGl6ZSBlbXB0eSB2ZWN0b3IgZm9yIGVpZ2VudmFsdWVzCkVJR0VORlVOIDwtIE5VTEwgICAgICAgIyBpbml0aWFsaXplIE5VTEwgZm9yIGVpZ2VuZnVuY3Rpb25zIG1hdHJpeApJTkRFWCA8LSBjKCkKCmZvciAoaiBpbiAwOm5fZmluaXRlKSB7CiAgICBsYW1iZGFfaiA8LSBrYXBwYV4yICsgKGoqcGkvMileMgogICAgZV9qIDwtIHRhZHBvbGUuZWlnKGosZ3JhcGgpJHBoaSMvc3FydCgxMDAwKQogICAgRUlHRU5WQUwgPC0gYyhFSUdFTlZBTCwgbGFtYmRhX2opICAgICAgICAgIyBhcHBlbmQgc2NhbGFyIHRvIHZlY3RvcgogICAgRUlHRU5GVU4gPC0gY2JpbmQoRUlHRU5GVU4sIGVfaikgICAgICAgICAgICAjIGFwcGVuZCBjb2x1bW4gdG8gbWF0cml4CiAgICBJTkRFWCA8LSBjKElOREVYLCBqKQogICAgaWYgKGo+MCAmJiAoaiAlJSAyID09IDApKSB7CiAgICAgIGxhbWJkYV9qIDwtIGthcHBhXjIgKyAoaipwaS8yKV4yCiAgICAgIGVfaiA8LSB0YWRwb2xlLmVpZyhqLGdyYXBoKSRwc2kjL3NxcnQoMTAwMCkKICAgICAgRUlHRU5WQUwgPC0gYyhFSUdFTlZBTCwgbGFtYmRhX2opICAgICAgICAgIyBhcHBlbmQgc2NhbGFyIHRvIHZlY3RvcgogICAgICBFSUdFTkZVTiA8LSBjYmluZChFSUdFTkZVTiwgZV9qKSAgICAgICAgICAgICMgYXBwZW5kIGNvbHVtbiB0byBtYXRyaXgKICAgICAgSU5ERVggPC0gYyhJTkRFWCwgaiswLjEpCiAgICB9Cn0KCiMgQ29tcHV0ZSBpbm5lciBwcm9kdWN0IG1hdHJpeApJTk5FUl9QUk9EVUNUUyA8LSB0KEVJR0VORlVOKSAlKiUgKEVJR0VORlVOICogd2VpZ2h0cykKCiMgQXNzdW1pbmcgSU5ORVJfUFJPRFVDVFMgaXMgYSBudW1lcmljIG1hdHJpeApwbG90X2x5KAogIHogPSBJTk5FUl9QUk9EVUNUUywKICB0eXBlID0gImhlYXRtYXAiLAogIGNvbG9yc2NhbGUgPSAiVmlyaWRpcyIgICMgWW91IGNhbiB0cnkgIkhvdCIsICJKZXQiLCBldGMuCikgJT4lCiAgbGF5b3V0KAogICAgdGl0bGUgPSAiSW5uZXIgUHJvZHVjdCBNYXRyaXggb2YgRWlnZW5mdW5jdGlvbnMiLAogICAgeGF4aXMgPSBsaXN0KHRpdGxlID0gIkluZGV4IiksCiAgICB5YXhpcyA9IGxpc3QodGl0bGUgPSAiSW5kZXgiKQogICkKCmBgYAoKCgpgYGB7cn0KRUlHRU5GVU4gPC0gYXBwbHkoRUlHRU5GVU4sIDIsIG9yZGVyX3RvX3Bsb3QpCnggPC0gb3JkZXJfdG9fcGxvdCh4KQp5IDwtIG9yZGVyX3RvX3Bsb3QoeSkKCiMgQ3JlYXRlIGludGVyYWN0aXZlIHBsb3QKZmlnIDwtIHBsb3RfbHkoeCA9IH5JTkRFWCwgeSA9IH5FSUdFTlZBTCwgdHlwZSA9ICdzY2F0dGVyJywgbW9kZSA9ICdtYXJrZXJzJywKICAgICAgICAgICAgICAgbWFya2VyID0gbGlzdChzaXplID0gNiksCiAgICAgICAgICAgICAgIG5hbWUgPSAiRWlnZW52YWx1ZXMiKQoKZmlnIDwtIGZpZyAlPiUgbGF5b3V0KHRpdGxlID0gIkVpZ2VudmFsdWVzIiwKICAgICAgICAgICAgICAgICAgICAgIHhheGlzID0gbGlzdCh0aXRsZSA9ICJpbmRleCBpIiksCiAgICAgICAgICAgICAgICAgICAgICB5YXhpcyA9IGxpc3QodGl0bGUgPSAiRWlnZW52YWx1ZXNfaSIpKQoKCgpwbG90X2RhdGEgPC0gZGF0YS5mcmFtZSgKICB4ID0gcmVwKHgsIHRpbWVzID0gbmNvbChFSUdFTkZVTikpLAogIHkgPSByZXAoeSwgdGltZXMgPSBuY29sKEVJR0VORlVOKSksCiAgeiA9IGFzLnZlY3RvcihFSUdFTkZVTiksCiAgZnJhbWUgPSByZXAoSU5ERVgsIGVhY2ggPSBsZW5ndGgoeCkpCikKCiMgQ29tcHV0ZSBheGlzIGxpbWl0cwp4X3JhbmdlIDwtIHJhbmdlKHgpCnlfcmFuZ2UgPC0gcmFuZ2UoeSkKel9yYW5nZSA8LSByYW5nZShFSUdFTkZVTikKCiMgSW5pdGlhbCBwbG90IHNldHVwIChmaXJzdCBmcmFtZSBvbmx5KQpwIDwtIHBsb3RfbHkocGxvdF9kYXRhLCBmcmFtZSA9IH5mcmFtZSkgJT4lCiAgYWRkX3RyYWNlKAogICAgeCA9IH54LCB5ID0gfnksIHogPSB+eiwKICAgIHR5cGUgPSAic2NhdHRlcjNkIiwgbW9kZSA9ICJsaW5lcyIsCiAgICBuYW1lID0gcGFzdGUwKCJuX2Zpbml0ZSA9ICIsIG5fZmluaXRlKSwKICAgIGxpbmUgPSBsaXN0KGNvbG9yID0gImJsdWUiLCB3aWR0aCA9IDIpCiAgKSAlPiUKICBsYXlvdXQoCiAgICBzY2VuZSA9IGxpc3QoCiAgICAgIHhheGlzID0gbGlzdCh0aXRsZSA9ICJ4IiwgcmFuZ2UgPSB4X3JhbmdlKSwKICAgICAgeWF4aXMgPSBsaXN0KHRpdGxlID0gInkiLCByYW5nZSA9IHlfcmFuZ2UpLAogICAgICB6YXhpcyA9IGxpc3QodGl0bGUgPSAiVmFsdWUiLCByYW5nZSA9IHpfcmFuZ2UpCiAgICApLAogICAgdXBkYXRlbWVudXMgPSBsaXN0KAogICAgICBsaXN0KAogICAgICAgIHR5cGUgPSAiYnV0dG9ucyIsIHNob3dhY3RpdmUgPSBGQUxTRSwKICAgICAgICBidXR0b25zID0gbGlzdCgKICAgICAgICAgIGxpc3QobGFiZWwgPSAiUGxheSIsIG1ldGhvZCA9ICJhbmltYXRlIiwKICAgICAgICAgICAgICAgYXJncyA9IGxpc3QoTlVMTCwgbGlzdChmcmFtZSA9IGxpc3QoZHVyYXRpb24gPSAxMDAsIHJlZHJhdyA9IFRSVUUpLCBmcm9tY3VycmVudCA9IFRSVUUpKSksCiAgICAgICAgICBsaXN0KGxhYmVsID0gIlBhdXNlIiwgbWV0aG9kID0gImFuaW1hdGUiLAogICAgICAgICAgICAgICBhcmdzID0gbGlzdChOVUxMLCBsaXN0KG1vZGUgPSAiaW1tZWRpYXRlIiwgZnJhbWUgPSBsaXN0KGR1cmF0aW9uID0gMCksIHJlZHJhdyA9IEZBTFNFKSkpCiAgICAgICAgKQogICAgICApCiAgICApLAogICAgdGl0bGUgPSAiVGltZTogMCIKICApCgojIENvbnZlcnQgdG8gcGxvdGx5IG9iamVjdCB3aXRoIGZyYW1lIGluZm8KcGIgPC0gcGxvdGx5X2J1aWxkKHApCgojIEluamVjdCBjdXN0b20gdGl0bGVzIGludG8gZWFjaCBmcmFtZQpmb3IgKGkgaW4gc2VxX2Fsb25nKHBiJHgkZnJhbWVzKSkgewogIHQgPC0gSU5ERVhbaV0KICBlcnIgPC0gc2lnbmlmKEVJR0VOVkFMW2ldLCA0KQogIHBiJHgkZnJhbWVzW1tpXV0kbGF5b3V0IDwtIGxpc3QodGl0bGUgPSBwYXN0ZTAoImk6ICIsIHQsICIgfCBFaWdlbnZhbHVlX2k6ICIsIGVycikpCn0KCgpmaWcgICMgRGlzcGxheSB0aGUgcGxvdApwYgoKYGBgCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoKCgoK